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We present results of three dimensional numerical simulations of the merger of unequal-mass binary 
neutron stars in full general relativity. A T-law equation of state P = (V — l)pe is adopted, where P, 
p, e, and V are the pressure, rest mass density, specific internal energy, and the adiabatic constant, 
respectively. We take F = 2 and the baryon rest-mass ratio Qm to be in the range 0.85-1. The 
typical grid size is (633, 633, 317) for (x, y, z) . We improve several implementations since the latest 
work. In the present code, the radiation reaction of gravitational waves is taken into account with a 
good accuracy. This fact enables us to follow the coalescence all the way from the late inspiral phase 
through the merger phase for which the transition is triggered by the radiation reaction. It is found 
that if the total rest-mass of the system is more than ~ 1.7 times of the maximum allowed rest-mass 
of spherical neutron stars, a black hole is formed after the merger irrespective of the mass ratios. 
The gravitational waveforms and outcomes in the merger of unequal-mass binaries are compared 
with those in equal-mass binaries. It is found that the disk mass around the so formed black holes 
increases with decreasing rest-mass ratios and decreases with increasing compactness of neutron 
stars. The merger process and the gravitational waveforms also depend strongly on the rest-mass 
ratios even for the range Qm = 0.85-1. 

I. INTRODUCTION 

Binary neutron stars such as the Hulse- Taylor binary pulsar [1] adiabatically inspiral as a result of the radiation 
reaction of gravitational waves, and eventually merge. In the most optimistic scenario, the latest statistical study 
suggests that such mergers may occur approximately once per year within a distance of about 30 Mpc [2]. Even the 
most conservative scenario predicts an event rate approximately one per year within a distance of about 400 Mpc 
[2]. This implies that the merger of binary neutron stars is one of the promising sources for kilometer-size laser 
interferometric detectors, such as LIGO, TAMA, GEO600, and VIRGO [3,4]. 

Interest has also been stimulated by a hypothesis about the central engine of 7-ray bursts (GRBs) [5]. Recently, it 
has been found that many GRBs are of cosmological origin [5]. In cosmological GRBs, the central sources must supply 
a large amount of the energy <; 10 51 ergs in a very short time scale (order of milliseconds to minutes) . Most GRB 
models involve a stellar system resulting in a stellar-mass rotating black hole and a massive disk of mass ~ 0. 1-1 Mq, 
which could supply a large amount of energy by neutrino processes or by extracting the rotational energy of the black 
hole. GRBs may be classified into two classes. One is a long burst for which the duration of the bursts is longer than 
~ 1 sec and typically ~ 10 sec, and the other is a short burst for which the duration is typically ~ 100 msec. It has 
been recently suggested that the merger of binary neutron stars is a possible progenitor to producing short bursts. 

Hydrodynamic simulations employing full general relativity provide the best approach for studying the merger of 
binary neutron stars. Over the last few years, numerical methods for solving coupled equations of the Einstein and 
hydrodynamic equations have been developed [6-10] and now such simulations are feasible. In previous papers [7,8], 
we focused on the binary neutron stars of equal mass, and have found the following results: (i) the final outcome (of 
either a neutron star or a black hole) depends on the compactness of each neutron star and on the equation of state. 
Even if the total mass of the system is ~ 1.5 times larger than the maximum allowed rest-mass of a spherical star 
for a given equation of state, a differentially rotating neutron star supported by a significant centrifugal force may 
be formed; (ii) in the case of neutron star formation, nonsphcrical oscillation modes of the formed neutron star are 
excited and, as a result, gravitational waves with characteristic frequency ~ 2-3 kHz are emitted; (iii) in the case of 
black hole formation, the disk mass around the formed black hole is negligible because the specific angular momentum 
of all the mass elements in equal-mass binary neutron stars is too small and also because the angular momentum 
transfer is not effective during the merger. 

So far, all the simulations in general relativity have been performed assuming that two neutron stars are identical 
[7,8,11,12], since they are indeed approximately identical in the observed systems of binary neutron stars [13]. For 
example, mass ratio of the Hulse- Taylor binary is about 0.963 [14]. However, it seems there is no theoretical reason 
that nature should produce only binary neutron stars of nearly equal mass. Allowed mass range of neutron stars may 
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fall in a fairly broad range ~ 1-2M Q according to theories of neutron stars [15,16]. From the theoretical point of 
view, it is reasonable and an interesting subject to investigate the merger of two unequal-mass neutron stars. 

One of the most important findings in the previous works [7,8] is that black hole formation is not accompanied 
by disks with a large mass. The mass of the disk is found to be less than 0.01M Q . This result suggests that binary 
neutron stars of equal mass may not be good progenitors for the central engine of GRBs. On the other hand, the disk 
mass may be much larger in the merger of binary neutron stars of unequal mass, because the smaller-mass neutron 
star would be tidally disrupted by the more massive primary before contact and would subsequently form a tidal tail 
around the primary in which angular momentum transfer is likely to be efficient to form disks around the central 
object. In addition, in association with the change of the merger process, gravitational waveforms may be significantly 
modified. Actually, Newtonian and post-Newtonian simulations indicate such significant changes [17-20]. 

From a computational point of view, we have substantially improved our implementation for a solution of Einstein 
and hydrodynamic equations from our previous approach [7,8]. Primarily, a modified numerical scheme for solving 
hydrodynamic equations by adopting the so-called high-resolution shock-capturing scheme [10,21] provides better 
accuracy. The spatial gauge condition is changed from a minimal distortion type [22,23] to a dynamical one in 
which a hyperbolic type equation is adopted for determining the shift vector [24,25]. This has resulted in saving 
a substantial amount of computation time. Finally, we have modified the treatment for the transport terms in the 
evolution equations of geometric variables. This improves the accuracies for a solution of the geometric quantities 
and for conservation of the total ADM mass and angular momentum significantly. 

The paper is organized as follows. In Sec. II, we review basic equations, gauge conditions, and methods for 
setting initial conditions currently adopted in fully general relativistic simulations of binary neutron star mergers. 
In Sec. Ill, methods used for analysis of gravitational waves are summarized. In Sec. IV, the numerical results 
are presented, paying particular attention to merger process, disk mass, and gravitational waveforms. Section V is 
devoted to a summary. Throughout this paper, we adopt the geometrical units in which G = c = 1 where G and c 
are the gravitational constant and the speed of light. Latin and Greek indices denote spatial components (x, y, z) and 
space-time components (t,x,y,z), respectively. Sij(— 5 l i) denotes the Kronecker delta. 



II. FORMULATION 
A. Basic equations 

In our numerical simulation, the Einstein and general relativistic hydrodynamic equations are solved without any 
approximation. The formulation for a numerical solution of these coupled equations is based on those described in 
a previous work [8]. However, we have improved several numerical implementations since then, and the form of the 
basic equations adopted in numerical simulation has also been modified. A summary of the current formulation is, 
therefore, in order here. 

The line element is written in the form 

ds 2 = g^dx^dx" = (-a 2 + f3 k f3 k )dt 2 + 2(3 l dx l dt + j^dx'dx 3 , (2.1) 

where g^, a, (3 l {fii = 'Yijft), and 7^- are the four dimensional spacetime metric, the lapse function, the shift vector, 
and the three dimensional spatial metric, respectively. Following [26,23,6], we define the quantities as 

7 = det( 7lJ )^e 1 ^, (2.2) 
~ ll3 =e-^ ll31 (2.3) 

M 3 ; = e-**{K ij -'^ 1ij K), (2.4) 

where Kij is the extrinsic curvature, and K its trace. In the Cartesian coordinates adopted in our simulation, det(7y) 
should be unity. In the numerical computations, <fi, jij, K, and A i3 are evolved in time, instead of 7^ and K^. 
Note that the indices of Aij (A'- 7 ') are raised (lowered) in terms of 7 y (7y). Hereafter, Di and Di are used as 
covariant derivatives with respect to 7^ and 7^, respectively. In addition, the Laplacians are defined as A = D l Di 
and A = D l Di. 

As the matter source of the Einstein equation, a perfect fluid is adopted. Then, the energy-momentum tensor is 
written as 

T„ v = (p + pe + P)u„u v + p 9^, (2-5) 
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where p, e, P, and u M are the baryon rest-mass density, the specific internal energy density, the pressure, and the 
four- velocity, respectively. Initial conditions are given using a polytropic equation of state as 



P = K P l 



r = i + - 

n 



(2.6) 



where k and n are a polytropic constant and a polytropic index. During the time evolution, we adopt a F-law equation 
of state of the form 



P = (T - l)pe. 



(2.7) 



In the absence of shocks, the polytropic form of the equation of state is preserved even if Eq. (2.7) is used. Thus 
the quantity «/ = Pj p r [= e/(T — l)p r_1 ] measures the efficiency of the shock heating. In this paper, we set n = 1 
(r = 2) as a reasonable qualitative approximation to moderately stiff equations of state for neutron stars [15]. 
The hydrodynamic equations (continuity, Euler, and energy equations) are written in the forms 



dtp* + d l {p*v % ) = 0, 
d t {p*Uj) + di(p*v l Uj + Pae 6 
d t {p*e) + d t [p*ev l + Pe 6 V 



? j ) = Pd j (ae**) 
-/?)] = ae^PK- 



p* 
P* 



whdjCt — UidjP 1 



1 



i*h 



3 llJ3r - ' 2u t h 
UiUjK lJ — p*Uij tJ Djct, 



th u k u l d j j kl 



where 



p* = pwe' v , 

P 

h = l + e+—, 
P 

w = au 1 , 
ii k = hu k , 

e 60 p 

e = Taun' J 'n u = hw , 

p* pw 



OTj %3 Uj 

whe 4 ^ ' 



(2.8) 
(2.9) 
(2.10) 

(2.11) 
(2.12) 

(2.13) 
(2.14) 

(2.15) 
(2.16) 



Here, the conservative form of the energy equation is adopted in contrast with the previous works [7,8]. In numerical 
simulations, Eqs. (2.8)-(2.10) are solved to evolve p*, and e. Once Uj is obtained, w is determined from the 



normalization relation of the four-velocity, u^u^ 



-1, which can be written as 



w 



1 + f lj UiUj = 1 + "/ tj UiU 3 



pw' 



(2.17) 



where P and p are related to p*, e, and w as P = P(p,e) = P[p*/(we 6 ^), e] and p = p*/(we 6 ^). 

The Einstein equation is split into the constraint and evolution equations. The Hamiltonian and momentum 
constraint equations are written in the form 



AV> = ^R k - 2tt Ph ^ 5 



(2.18) 
(2.19) 



where tp = e^, pn = T liV n^n Vl and J; = —T^n^vi with = (— a,0). Here, Rij (Rij) denotes the Ricci tensor 
with respect to jij (7^), and R k k = Rijj^ (R k k = Rijj^). These constraint equations are solved only at t — to set 
initial conditions (see Sec. II D) and for t > 0, they are used to monitor the accuracy of numerical solutions. 
Following [26,23,6-8], evolution equations for the geometric variables are written as 



(2.20) 



(d t -p l d l )A ij =e 



- o -40 
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+a{KA ij - 2A ik Af) + p\ t A kj + j3 k /A kl - ^ k A tJ 

^a(e-^S l3 - l -% 3 S k k ), (2.21) 
(fit - fidfiK - a[A l0 A^ + i^ 2 ] - Aa + 47ra(p H + S k k ), (2.22) 

where Sij = T^j^-y^j. Equations (2.20)-(2.22) are solved to evolve -yy, Aij, and K. 
In the previous works [26,23,6-8], the evolution equation for <j> is written in the form 

{d t -p l d l )^= l -(-aK + (3 k tk ). (2.23) 

Instead of this form, in the present work, a conservative form is adopted as 

8 t e^ - d % {f3 l e^) = -uKe H . (2.24) 

As in the previous works, we introduce an auxiliary variable Fi = S^dfiij [26], which evolves according to the 
evolution equation 

{dt - tfdfiFi = -WiraJi + 2a{f kj A ikJ + f k ^A lk - \ti l h jhi + 6</> fe i* - ^K ti } 

+ ^{-20,^+^^,1 + (pratfj - \lu&,!) J, (2-25) 

where jij and 7 y are split into Sij + hij and (5 y + / u . In the numerical simulations, a term 5 kl ji k ,ij which appears 
in the expression of Rij in Eq. (2.21) is evaluated using Fi as Fij. This replacement is crucial to enable a stable and 
longterm simulation. 

B. Improvements for numerical implementations 

In this section, we report our improvement to several numerical implementations since the latest work [8] . Firstly, the 
hydrodynamic part has been improved adopting the so-called high-resolution shock-capturing scheme for computation 
of the transport terms of Eqs. (2.8)-(2.10) as described in [21]. With this scheme, shocks are captured with a much 
better accuracy than in the previous implementation [6-8]. Although the shocks generated during the merger of 
binary neutron stars are not very strong, it is promising to use such high-resolution schemes to accurately compute 
peak densities and to evaluate effects of the shock heating. 

Numerical treatments for transport terms in the evolution equation for geometric variables have been also improved. 
For the transport terms in Eqs. (2.20)-(2.22), a second-order upwind scheme has been adopted [23]. To avoid numerical 
instabilities, we incorporate a limiter / by which the order of the finite differencing for the numerical flux is lowered 
from the second order to the first order at a point of steep gradient as 

F mm ^F 1 f + F 2 (l-f), (2.26) 

where F num , F\, and F 2 denote the total, first-order, and second-order fluxes. Since the previous choice of / is found 
to be too dissipative [23], we have changed the functional form of / to 



2y/\SQJQ d \ 

J \SQ u \ + \SQ d \ 7 [ • ' 

where Q denotes one of the variables among A\j, K, and 7^ . 5Q U and SQd denote the difference of Q for two 
neighboring grid points as SQ U = Q1+1 — Qi and SQd = Qj — Q7-1 where Qi denotes the value of Q at the 7-th grid 
point. 

For the evolution of e 6 ^, we have also changed the finite differencing scheme. As shown in Eq. (2.24), the evolution 
equation for e 6 ^ has the same conservative form as that of hydrodynamic equations. Thus the numerical flux is 
computed using the third-order upwind scheme with an appropriate min-mod limiter as done in the hydrodynamic 
equations [21]. This change plays a significant role for enforcing the conservation of the total ADM mass and angular 
momentum. 
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The outer boundary condition for <j> has been also improved. In our previous works, we simply imposed 

< j> = 0(r- 1 ). (2.28) 

It is replaced with a better condition as 

<^=^ + 0(r- 2 ), (2.29) 
where M is the ADM mass which will be defined in Sec. II E. 

C. Change in gauge conditions 

As the time slicing condition, an approximate maximal slice (AMS) condition K w is adopted following previous 
papers [8]. On the other hand, the spatial gauge condition has been changed. 

Our previous simulations were performed adopting an approximately minimal distortion (AMD) gauge condition 
[23] . The equation in this condition is schematically written as 

5 tf A f /J* + hdidtP = S u (2.30) 

where Si denotes the source term composed of geometric variables and matter sources, and A/ the flat Laplacian. 
In this gauge condition, a vector elliptic-type equation has to be solved. A serious drawback of this is the long 
computational time needed to obtain its numerical solution. Typically, <~ 50% of total computational time is consumed 
in solving this equation. 

To overcome this drawback, we adopt a dynamical spatial gauge condition, e.g., in [24,25]. Following [27], the 
equation for the shift vector is chosen to be 

d t (i k = j kl [F l +At(d t F l )], (2.31) 

where At denotes a time step in numerical computation. The second term in the right-hand side of Eq. (2.31) 
is introduced to stabilize the numerical computation. In this choice, (3 k obeys a hyperbolic type equation (for a 
sufficiently small value of At) as 

d 2 t F = A^ + ^dkdiF - S j , (2.32) 

where S J denotes the source term. With this gauge condition, the fraction of the computational time occupied for 
imposing the spatial gauge is negligible. Furthermore, we have confirmed that the numerical solution in this gauge 
condition agrees well with that in the AMD gauge condition for collapse of neutron stars to a black hole [27] and 
for oscillating and rapidly rotating neutron stars (cf. Fig. 1). This indicates that the present dynamical gauge is 
physically as good as the AMD gauge. 

Since (3 k obeys a hyperbolic type equation in this gauge condition, so should F;. Thus we impose an outgoing 
boundary condition for Fi as 

Fl = W- r \ (2.33) 
r 

where /; (t — r) is a function and its value at outer boundaries is determined from the values on the eight nearby grid 
points at a previous time step. The same type of boundary condition is imposed for and 7^ [26] . 
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FIG. 1. Evolution of central density p c in units of the initial value p c ,o and central value of the lapse function a c for an 
oscillating and rapidly rotating neutron star with n = 1. The baryon rest-mass and the ADM mass in units of n = 1 (see 
Sec. II E) are 0.186 and 0.172, respectively. The compactness measured by the equatorial (polar) radius is 0.129 (0.207). The 
angular velocity is constant and equal to the Kepler velocity at the equatorial surface. The time is shown in units of p c . 
The solid and dashed curves denote the results by the dynamical and AMD gauge conditions, respectively. 

D. Initial conditions 

Binary neutron stars with a moderate compactness of orbits as a/M ^ 6 where a denotes an orbital separation 
are in a quasicquilibrium state even just before the merger because the time scale of gravitational radiation reaction 
at Newtonian order ~ 5/{64£l(MNf2) 5//3 } [15] (where Mn and Vl denote the Newtonian total mass of system and the 
orbital angular velocity of binary neutron stars) is several times longer than the orbital period. Thus a quasiequilibrium 
state should be prepared as the initial condition for a realistic simulation of the merger. Such quasicquilibrium states 
are obtained by solving coupled equations of gravitational field and hydrostatic equations. For the gravitational 
field, we adopted the conformal flatness formulation [28] in which the three geometry is assumed to be conformally 
flat and the selected components of the Einstein equation are solved. Specifically, the selected components are the 
Hamiltonian and momentum constraints, and the trace of spatial projection of the Einstein equation with the maximal 
slicing condition K = 0. The solution in this formalism is fully general relativistic in the sense that they satisfy the 
constraints. 

It is expected that most of the close binary neutron stars in quasiequilibrium circular orbits have irrotational velocity 
fields approximately since the viscous time scale is much longer than the gravitational radiation time scale and the 
orbital period ~ 2 msec is much shorter than the typical spin period of neutron stars [29] . Assuming the irrotational 
velocity field and the presence of a helical Killing vector as 




(2.34) 



the hydrodynamic equations are written into a first integral of the Euler equation and an elliptic-type equation for a 
velocity potential [30]. 

The coupled equations of the selected Einstein and hydrostatic equations are solved by a pseudospectral method 
developed by Bonazzola, Gourgoulhon, and Marck [31]. Detailed numerical calculations have been done by Taniguchi 
and part of the numerical results are presented in [32]. 

Quasiequilibrium solutions are given as the initial conditions for simulations without any modification. In a realistic 
system of binary neutron stars, the orbit is not strictly circular because of the presence of the approaching velocity 
due to gravitational radiation reaction. As pointed out by Miller [33], neglecting the effect of the approaching velocity 
yields a systematic error in waveforms and the merger process, if we choose a quasiequilibrium with a very small 
orbital separation as the initial condition. It is likely that gravitational waveforms, in particular the wave phase, 
obtained below contain a small systematic error. However, it has been studied in [8] that the merger process and the 
final outcome depend very weakly on an artificial approaching velocity of sa 10% of the orbital velocity. 
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E. Definitions of quantities 



In numerical simulations, we refer to the total baryon rest-mass, the ADM mass, and the angular momentum of 
the system, which are given by 



p*d 3 x, 



M=-—<t D l iPdSi 
2tt 



Rue 



,5d> 



16tt 



AijA ij - -K 2 - R k k e~ 44> 



(fx, 



J= 



1 

8^ 



d 6 x, 



(2.35) 



(2.36) 



(2.37) 



where dSj = r 2 Djrd(cos9)dip and ^ — —y(d x y + x(d y y . To rewrite the expressions for M and J, the Gauss law is 
used. Here, M* is a conserved quantity, and it uniquely specifies a model of a stable neutron star for a given value of 

r. 

M and J are computed using the volume integral shown in Eqs. (2.36) and (2.37). Since the computational domain 
is finite, they are not constant and decrease after gravitational waves propagate to the outside of the computational 
domain during time evolution. Therefore, in the following, they are referred to as the ADM mass and the angular 
momentum computed in the finite domain (or simply as M and J, which decrease with time). As easily predicted from 
the calculation using the quadrupole formula, M decreases at most by 0.5% and may be regarded as an approximately 
conserved quantity, while J decreases by ~ 5-10%. 

A model of each neutron star is specified using the compactness (M/R)^ which is defined as the ratio of the ADM 
mass to the circumferential radius of a spherical neutron star in isolation (see Tables I and II) . To indicate how massive 
the system is, we also introduce the ratio of the total baryon rest-mass of the system to the maximum allowed mass 
of spherical neutron stars for a given equation of state M*max, 



Q* = 



M sph 



(2.38) 



Physical units enter the problem through the polytropic constant k initially chosen, which can be completely scaled 
out of the problem. Since n 11 / 2 has the dimension of length, time, and mass in the geometrical units c = G = 1, 
dimcnsionless variables can be constructed as 



M* = M*k - "/ 2 , M = MK~ n/2 , R = R K - n/2 , 
J=jK~ n , p = pn n , and ti = n K n / 2 . 



(2.39) 



In the following, only these dimensionless quantities are presented (namely units of k = 1 are adopted) and, hence, 
the bar is omitted. 

Nondimcnsional quantities may be converted to dimensional ones for a value of k. For k = 1.455 x 10 cgs which 
is chosen in [9], the mass, the density, and time in the dimensional units are written as 



M dim = 1.80M C 



1/2 



M 



1.455 x 10 5 cgs/ V - 180 / 

-l 

Pdim - 1.86 x 10 15 g/cm 3 ' 



1.455 x 10 5 cgs/ V - 300 / 



Tdim = 4.93 msec 



1.455 x 10 5 cgs 



1/2 



T 
100 



(2.40) 
(2.41) 
(2.42) 
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F. Calibration 



Several test simulations, including spherical collapse of dust, stability of spherical and rotating neutron stars, 
comparison of eigenoscillation modes of spherical stars with the known results, and longterm evolution of rotating 
stars, have been performed to check the reliability of numerical results obtained in the new implementation. A list of 
these test simulations and some of their results are described in [21]. 

During the simulations, we monitored the violation of the Hamiltonian constraint, and the conservation of the 
baryon rest-mass, the ADM mass, and the angular momentum. Because of the emission of gravitational waves, M 
and J computed in the finite volume by Eqs. (2.36) and (2.37) decrease with time. However, the sum of M and 
accumulated radiated energy of gravitational waves, and the sum of J and accumulated radiated angular momentum 
of gravitational waves should be conserved (at least approximately) in numerical computation as [34] 

M(t) + AE(t) = M , (2.43) 
J(t) + AJ{t) = J , (2.44) 

where AE(t) and A J(t) denote the total radiated energy and angular momentum by gravitational waves until a time 
t, for which the definitions are described in Sec. III. Mo and Jo denote the initial values of M and J. 
The violation of the Hamiltonian constraint is locally measured by the equation as 



Ay, - ±k k k + 2 w 5 + g (jgjft - fif 2 ) 



(2.45) 



IAVI + If i? fe fe | + |2W 5 I + %(\A ij AV\ + IK*) ' 
Following [21], we define and monitor a global quantity as 

H=^-J P*Ud 3 x. (2.46) 
Hereafter, this quantity will be referred to as the averaged violation of the Hamiltonian constraint. 

III. ANALYSIS OF GRAVITATIONAL WAVES 

Gravitational waves are measured in terms of the gauge-invariant Moncrief variables in a flat spacetime [35]. To 
compute them, first, we perform a coordinate transformation for the three-metric from the Cartesian coordinates to 
the spherical polar coordinates, and then split 7^ into rjij + X^mCiJ 1 ; where r]ij is the flat metric in the spherical 
polar coordinates and (Ij 1 is given by 

r 2 (K lm Y lm + G lm W lm ) r 2 G lm X lm 
* * r 2 sin 2 9(K lm Y lm - Gi m Wi m ) 

-Ci m d v Yi m / sm6 Ci m d Yi m sm8 \ 

* r 2 D lm X lm /sm8 -r 2 D lm W lm sm6 . (3.1) 

* * -r 2 D lm X lm sin 6 J 

Here, * denotes the symmetric components. The quantities H% m , hu m , Ki m , Gi m , Ci m , and D; m are functions of r 
and t, and are calculated by performing integrals over a two-sphere of a given coordinate radius [see [26] for details]. 
Y\. m is the spherical harmonic function, and Wi m and Xi m are defined as 



W lr 



(d e ) 2 - cotOde - —^(d v ) 2 Y lm , X lm = 2d v d e - cote Y lm . (3.2) 



1 

sh?6> 

The gauge-invariant variables of even and odd parities are defined by 



R? m (t,r) = \Jjf^jr{^2i m + l(l + l)k Um }, (3.3) 
R? m (t,r) = J^3( ^+rd r D lm ), (3.4) 



(1-2)! V r 
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where 



k Um = K lm + 1(1 + l)Gi m + 2rd r Gi m - 2 



hlln 



H. 



lira 



id_ 

2dr 



r{K lm + 1(1 + l)G im } 



The cosine and sine components of the gauge-invariant variables, which are real quantities, are also defined as 



R E + i? E 
Ri m + = — and Rirr, 



pE _ pE 
"-Im n l -m 



V2 



V2 



(m > 0). 



(3.5) 
(3.6) 

(3.7) 



Using the gauge-invariant variables, the energy luminosity and the angular momentum flux of gravitational waves 
can be calculated as 



f-a F E[i*«ti' + iw«. 

l,m 



Im 



l.ra 



The total radiated energy and angular momentum are defined as 



(3.8) 
(3.9) 

(3.10) 



We have computed modes with I = 2, 3, and 4, and found that the even modes with I — \m\ — 2 are dominant, and 
the even mode with I = 2 and m = is secondly dominant. For merger of unequal-mass binaries, the amplitude of the 
even modes with I = \m\ = 3 are as large as that of I = 2 and m = . Thus attention is paid to these three modes. 
To search for the dominant frequencies of gravitational waves, the Fourier spectra are computed by 



Ri 



(/, =I 



e 2mft R lm ±dt. 



(3.11) 



In the analysis, tf is chosen as the time at which the simulation is stopped. Before t < r Q b s where r b s denotes a 
radius at which gravitational waves are extracted, no waves propagate to r Q b s , so that we choose w r Q b s . 

Using the Fourier spectrum, the energy spectrum which is often referred to in literature (e.g., [18,20]) can be written 

as 



d 4 = YY. 

- 1 l,m>0 



where for m ^ 0, wc define 



Rim = ^\Rl m +(f)\ 2 + \Rl m -(f)\ 2 . 

To help the calculation of dE/df, \Rim(f)f\r is presented as the Fourier spectrum in the following. 



(3.12) 



(3.13) 



IV. NUMERICAL RESULTS 
A. Setup for simulation 

Several quantities that characterize quasiequilibrium states of irrotational binary neutron stars used as initial 
conditions for the present simulations are summarized in Table I. All quantities are appropriately scaled with respect 
to k to be dimensionless. 

As the initial conditions, we choose binaries of an orbital separation which is slightly (by ~ 10%) larger than that 
for an innermost orbit. Here, the innermost orbit is defined as a close orbit for which Lagrange points appear at the 
inner edge of neutron stars [36,31]. Models M1414 and M1616 are equal-mass binaries, and others are unequal-mass 
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ones. Total baryon rest-masses for models M1616, M1517, and M1418 or for M1414 and M1315 are almost identical, 
while the baryon rest-mass ratios are identical for models M1517 and M159183 as 0.925. 
The frequency of gravitational waves for binaries in these quasiequilibria is given by 

and, thus, the orbital period of the quasiequilibria, Pt=o, is 

/2.8M VY Co Y 3/2 , s 

P« « 2.08 msec _° ^ , (4.2) 



where Co is a compactness parameter of an orbit defined by 

Co = (M O) 2 /3 = Ml. (4.3) 
ao 

Here, ao is defined as the initial orbital separation. For the initial conditions chosen in this paper, ao > 8.5Mo. 

The simulations were performed using a fixed uniform grid and assuming reflection symmetry with respect to the 
equatorial plane (here, the equatorial plane is chosen as the orbital plane). The typical grid size is (633, 633, 317) 
for (x, y, z). The grid covers the region — L < x < L, —L < y < L, and < z < L where L is a constant. The grid 
spacing (which is L/316 in the typical case) is determined from the condition that the major diameter of each star is 
covered with about 40 grid points initially. 

Numerical results depend weakly on the grid resolution and location of the outer boundaries. In order to investigate 
this, additional test simulations were performed choosing the smaller grid sizes with a fixed value of grid spacing and 
the larger grid spacings with a fixed value of L for selected models. The setting for the test simulations are summarized 
in Table II and the numerical results are presented in Sec. IV E. 

With a (633, 633, 317) grid size, about 240 GBytes computational memory is required. For the case of neutron star 
formation, the simulations are performed for about 20000 time steps and then stopped artificially. The computational 
time for one model in such a calculation is about 100 CPU hours using 32 processors on FACOM VPP5000 in the 
data processing center of National Astronomical Observatory of Japan (NAOJ). For the case of black hole formation, 
the simulations crash soon after the formation of apparent horizon because of the so-called grid stretching around the 
black hole formation region. In this case, the computational time is about 50 CPU hours for about 10000 time steps. 

In the above setting, the wavelength of gravitational waves at t = (denoted by Ao) is about twice that of L (cf. 
Table I). As found in a previous paper [8], gravitational waves and radiation reaction are taken into account with 
a fair accuracy (within ~ 10% numerical error) in this setting. Since the typical wavelength of gravitational waves 
becomes shorter and shorter in the late inspiral phase, the accuracy of the wave extraction is improved with the 
evolution of the system. As a result, the magnitude of the error in the total radiated energy and angular momentum 
would be much smaller than 10%. This point will be reconfirmed in Sec. IV D. The wavelength of quasiperiodic waves 
emitted from the formed neutron star is much shorter than Ao and L, so that the waveforms in the merger stage can 
be computed accurately in the case of neutron star formation. 

As found in [36,31], orbits for irrotational binaries of equal mass with T < 2.5 (n > 2/3) are dynamically stable 
from the infinite separation to the innermost orbit. Therefore the merger in reality should be triggered by the 
radiation reaction of gravitational waves for T = 2. In the previous implementation, however, the radiation reaction 
for the late inspiral stage was not very accurately computed. Thus the simulations were initiated by setting a binary 
at the innermost orbit and reducing the angular momentum slightly to induce prompt merger [7,8]. In the new 
implementation, on the other hand, the radiation reaction can be computed with a good accuracy (within ~ 1-2 % 
error throughout a simulation, see Sec. IV D). Thus, in the present work, we prepared binaries of orbits slightly 
far away from the innermost orbits and started simulations without adding any perturbation. With this setting, a 
transition from the inspiral to the merger is triggered by the radiation reaction. This point will be demonstrated in 
Sec. IV D. 
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B. Characteristics of the merger 
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FIG. 2. Snapshots of the density contour curves for p in the equatorial plane for model M1315. The solid contour curves 
are drawn for p/0.15 = 1 — O.lj for j = 0, 1, 2, • • • , 9, and the dashed-solid curves are for p/0.15 = 0.05, 0.01, 10 -3 and 10 -4 . 
Vectors indicate the local velocity field (v x , v v ), and the scale is shown in the upper right-hand corner. Pt=o denotes the orbital 
period of the quasiequilibrium configuration given at t = 0. The length scale is shown in units of GMo/c 2 , where Mo is the 
gravitational mass computed at t — 0. In the first panel, the primary neutron star is located at x > 0. 
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FIG. 3. The same as Fig. 2 but for model M1414. 



In Figs. 2-5, we display the snapshots of the density contour curves and the velocity vectors at selected time steps 
for models M1315, M1414, M1418, and M1616 [37]. Figure 6 shows the evolution of the maximum values of p and </>, 
and the minimum value of a for all the models adopted in this paper except for model M159183. 

The numerical results for models M1414 and M1616 computed by an old implementation have been already presented 
in [7,8]. With the present new implementation, however, the quality of the numerical results is significantly improved 
and, hence, the improved results are displayed as the updated ones. 

The simulations for models M1616, M1517, and M1418 crashed soon after the formation of apparent horizons 
because the black hole forming region was stretched significantly and the grid resolution became too poor to resolve 
such a region. On the other hand, we artificially stopped the simulations for models M1414 and M1315 at t <; 3.5Pt = o 
to save the computational time. At the termination of these simulations, the averaged violation of the Hamiltonian 
constraint does not increase rapidly and remains of order 0.1. Therefore the simulations could be continued for a 
much longer time than 3.5P t= o in the formation of the massive neutron stars. 

In every model, the merger is triggered by the radiation reaction: For t < Pt=o, the orbital separation decreases 
gradually as a result of gravitational radiation reaction and each neutron star is elongated little by little. The 
elongation is always larger for the smaller-mass star in nonequal mass binaries. At a critical separation which is 
reached at t ~ 0.8-0.9Pt = o, the orbit becomes unstable probably against hydrodynamic instability to start the 
merger. At this point, the lag angle which is defined to be the angle in the equatorial plane between the major axis 
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of each star and the axis connecting the centers of mass of two stars [20] is ~ 10-15 degree. In unequal-mass binaries, 
we always find a larger lag angle for the smaller-mass star. This is in agreement with that in [20]. 

In contrast to previous Newtonian [17,18], post-Newtonian [19,20], and approximately rclativistic simulations [38], 
the formation of a black hole can be determined in fully general relativistic simulations. In the black hole formation, 
most of the fluid elements are swallowed into the black hole. Therefore the evolution of the system is significantly 
different from that in the neutron star formation. From this reason, we describe the character of the merger for the 
formation of neutron stars and black holes separately. 

1. Formation of hypermassive neutron star 

For models M1315 and M1414, massive neutron stars are formed (cf. Figs. 2 and 3). In this case, the total baryon 
rest-mass of the system is about 1.62 times as large as M*max for a polytropic equation of state with n = 1. Thus the 
formed neutron star is hypermassive in the sense that the mass is larger than the maximum allowed value for rigidly 
rotating neutron stars with n = 1 [39]. As indicated in the previous papers [7,8], such a large mass is supported by a 
large centrifugal force due to rapid and differential rotation. 

In [7,8], we concluded that the merged object for model M1414 eventually collapses to a black hole in 2-3 Pt=o- 
However, in the present improved simulation, a hypermassive neutron star is formed instead of a black hole. There 
are two plausible reasons for this discrepancy. One is that shocks are calculated in a better accuracy with the new 
implementation of a high-resolution shock-capturing scheme and, as a result, the thermal energy, which could play 
an important role for supporting the massive object, is increased in the new result. The other possible reason is an 
improvement on the treatment of the transport term for geometric variables. This makes the angular momentum 
conservation more accurate, avoiding the spurious collapse. 

Besides the correction to the threshold for collapse of merged objects, the qualitative properties during the merger 
of equal-mass binaries are essentially the same as those found in the previous simulations [7,8]: The merged object 
constitutes a double core structure soon after the merger [cf. fifth panel of Fig. 3 and Fig. 7(b)]. At the collision 
of two neutron stars, the radial infall velocity is not so large that shock heating is not very important around the 
mass center. In Fig. 8, we show k' along the x and y axes, which denotes the efficiency of the shock heating. In our 
notation, it is unity everywhere inside the neutron stars at t = 0. Thus the fluid elements for which the value of k' is 
larger than unity have experienced the shock heating. This figure shows that for r <J 4Mq, k' is ~ 1-2, implying that 
the shocks do not play a very important role except for the outer envelops. 

In the outer region, small spiral arms are formed soon after the merger sets in, but they do not spread outward 
widely because of insufficient angular momentum. The spiral arms subsequently wind around a central double core. In 
the central region, the double core has not only a rotational motion around the mass center but also has a quasiradial 
oscillation which is originally excited by a radial plunge at a transition from the inspiral to the merger stage. Because 
of the quasiradial oscillation, weak shocks are formed in the outer envelops. As a result, the outer region is heated up 
and gains the kinetic energy to expand outward [cf. Fig. 8(b)]. This process is repeated many times transferring the 
kinetic energy of the inner core to the outer region and, hence, the quasiradial oscillation of the core damps gradually 
(see Fig. 6). 

For the merger of an unequal-mass binary (for model M1315), the merger process is qualitatively different from 
that in the equal-mass case because tidal disruption of the smaller-mass star by the massive primary takes place (cf. 
fourth panel of Fig. 2). The tidally disrupted star subsequently forms a tidal tail. During the formation of the tidal 
tail, the angular momentum is efficiently transfered outward and, as a result, large spiral arms are formed. The spiral 
arms subsequently wind around the central core to be accretion disks. Because of the angular momentum transfer at 
the tidal disruption and at subsequent formation of spiral arms, the disk radius is much larger than that for model 
M1414 (cf. Fig. 9). 

In the central region, a massive object with an asymmetric double core is formed [cf. the last panel of Fig. 2 and 
Fig. 7(a)]. As in model M1414, the central core oscillates quasiradially (see Fig. 6). This motion produces shocks 
around the outer part of the core and, as a result, the energy is transfered to the outer envelops [cf. Fig. 8(a)]. Since 
this process is repeated, the quasiradial oscillation of the core damps gradually. The amplitude of the quasiradial 
oscillation for model M1315 is not as large as that for model M1414. This reflects the difference of the merger process 
between M1414 and M1315: For M1414, two neutron stars merge without tidal disruption and mass ejection outward. 
Therefore all the mass elements in this system collide coherently. On the other hand, for M1315, the tidal disruption 
takes place. As a result, a fraction of mass elements in the smaller-mass star do not have a plunging motion at 
the collision and, hence, the merger does not set in as coherent as that for model M1414. Due to this reason, the 
amplitude of the quasiradial oscillation is suppressed. 
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Figure 6(a) shows that the maximum density of the hypermassive neutron star for model M1315 is larger than 
that for M1414 in spite of the fact that the total baryon-rest mass is nearly identical. This reflects the fact that the 
region around the mass center for model M1315 rotates less rapidly than that for M1414. This suggests that the 
hypermassive neutron stars formed from the merger of the smaller rest-mass ratios are more compact. 

In Fig. 9, we plot the evolution of the baryon rest-mass fraction outside the spheres of radius 3M (solid curve), 
4.5Mo (dashed curve) and 6Mo (dotted-dashed curve). Here, r = is chosen as the center of the spheres. This shows 
the significance of the angular momentum transfer for model M1315. For model M1414, the baryon rest-mass outside 
the spheres of fixed radii simply oscillates with a mean value which is approximately constant with time evolution. 
The fraction of the rest-mass outside the sphere of radius 6Mq is <~ 1% in this case. On the other hand, for M1315, the 
baryon rest-mass outside spheres of fixed radii increases gradually with time. This result reflects an efficient angular 
momentum transfer. Figure 9 indicates that the fraction of the rest-mass outside the sphere of radius 6M is ~ 5%, 
implying that disks of J> O.IMq are formed around the hypermassive neutron star. 

A post-Newtonian simulation reports that the fraction of the disk mass for an equal-mass merger is ~ 6% [20] . 
This value is much larger than the value obtained in this paper. A plausible reason for this discrepancy is that in the 
post-Newtonian approximation, the gravity of the hypermassive neutron star is underestimated and, hence, the mass 
captured by it is also underestimated. 

Because of the nonaxisymmetric and quasiradial oscillations of the hypermassive neutron stars, quasiperiodic grav- 
itational waves of a few characteristic oscillation modes are simultaneously excited for models M1315 and M1414 for 
a long duration after the merger. This point will be discussed in Sec. IV C. 

2. Formation of rotating black hole 

For models M1418, M1517, M1616, and M159183, black holes are formed (cf. Figs. 4, 5, and 10) in a dynamical time 
scale ~ 1.3P t= o irrespective of the baryon rest-mass ratios. The formation of the black holes is determined by finding 
the apparent horizons [40]. In all the cases, the total baryon rest-mass of the system is about 1.75 times as large as 
M*max for a polytropic equation of state with n = 1 . Since the black holes are formed for M* £ 1.75M* s P h ax while the 
hypermassive neutron stars are the outcomes for <, 1.65M*£iax, the threshold of the total baryon rest-mass for 
the prompt black hole formation is between 1.65 and 1.75M 4 s it for n = 1. 

The formation process of the black holes depends on the rest-mass ratios. For the merger of two equal-mass neutron 
stars, the merger results in a massive object of a double core without tidal disruption and mass ejection outward. 
The merged object is too massive to support its self-gravity and, hence, collapses to a black hole promptly. Since the 
specific angular momentum of each fluid element is too small [8] and also since there is no efficient transfer of angular 
momentum during the merger, the disk mass around the formed black hole is very small [cf. Figs. 10(a) and (e)]. 

On the other hand, for the merger of two unequal-mass neutron stars, the black hole formation appears to be 
triggered by accretion to the primary star: First, the primary star tidally disrupts the smaller-mass companion at a 
critical separation. Subsequently, most of the tidal debris accrete to the massive primary star and a small fraction of 
them form spiral arms. The accretion increases the mass of the primary star rapidly, eventually, exceeding the critical 
value for formation of a black hole. During the merger, the angular momentum transfer works efficiently in the spiral 
arms, subsequently forming an accretion disk around the formed black hole. 
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FIG. 4. The same as Fig. 2 but for model M1418. Here, the solid contour curves are drawn for p/0.20 = 1 — O.lj for 
j = 0, 1, 2, • ■ ■ , 9, and the dashed-solid curves are for p/0.20 = 0.05, 0.01, 10~ 3 and 10" 4 . The thick dotted circle in the last 
panel of radius r ~ 0.5Mq denotes the location of the apparent horizon. At t — 0, the primary neutron star is located at x > 0. 
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(a) (b) 

FIG. 6. Evolution of the maximum values of p and <f), and the minimum value of a (a) for models M1414 (dashed curves) 
and M1315 (solid curves), and (b) for models M1616 (dashed curves), M1517 (long-dashed curves), and M1418 (solid curves). 
Note that for the case of black hole formation [Fig. 6(b)], the maximum density decreases in the final stage. The reason is as 
follows: We choose p* as a fundamental variable to be evolved and compute p from p»/iu/e 6< *. In the final stage, <j> is very large 
(> 1) and, hence, a small error in <f) results in a large error in p. Note that the maximum value of p* increases monotonically 
by many orders of magnitude. 
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FIG. 9. Evolution of the baryon rest-mass fraction outside the spheres of radius 3M (solid curve), 4.5Mo (dashed curve) 
and 6M (dotted-dashed curve) (a) for model M1315 and (b) for model M1414. 

In Fig. 10, we compare snapshots of density contour curves soon after the formation of apparent horizons for 
models M1616, M1517, M1418, and M159183. Obviously, fractions of the disk mass and the disk radius are larger for 
binaries of the smaller rest-mass ratios. Comparing the figures for M1517 and M159183 for which the rest-mass ratio 
is identical as 0.925, it is also found that the disk mass and radius are smaller for the merger of the larger compactness. 

Figure 11 shows the evolution of the baryon rest-mass fraction outside spheres of fixed coordinate radii for models 
M1616, M1517, M1418, and M159183. As the coordinate radii of the spheres, it is desirable to choose the radius of 
the innermost stable circular orbit around the formed black holes, but in practice it is difficult to determine it from 
numerical results exactly for dynamical spacetimes. Thus we rather arbitrarily choose 3M and 4.5M which are not 
far away from the radius of the innermost stable circular orbit for rotating black holes of nondimcnsional angular 
momentum parameter ~ 0.8-0.9. For model M1616, the mass fraction outside these radii decreases monotonically, 
and at the termination of the simulation, less than 0.2% of the total rest-mass of the system is outside the sphere 
of radius 3M . On the other hand, the mass fraction for r > 3M appears to approach ~ 2% and ~ 4% for models 
M1517 and M1418, respectively. This indicates that for mergers of unequal-mass neutron stars, a certain fraction of 
the mass would form disks, and the mass fraction seems to increase in proportion to 1 — Qm for a given value of M*. 

Comparing the results for models M1517 and M159183 for which the rest-mass ratio is identical as 0.925, it is 
found that the mass fraction of disks decreases with the increase of the compactness of neutron stars. The reason is 
simply that the gravity of the system is stronger for model M159183 and, as a result, a larger fraction of the mass is 
swallowed into the black hole. Obviously, smaller compactness of progenitor neutron stars is in favor of the formation 
of disks around a formed black hole. 

Figures 12(a)-(d) show n' along the x axis for models M1616, M1517, M1418, and M159183. Figures indicate that 
most of the fluid elements arc heated up by shocks, except the inner region of the disk where the value of k' is less 
than 10. This implies that the shock heating is not very important for high density regions; i.e., at the collision of 
two neutron stars, the shocks are not very strong. 
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FIG. 10. Comparison of the density contour curves for p in the equatorial plane near the end of the simulations for models 
M1616, M1517, M1418, and M159183. The solid contour curves and the velocity vectors are drawn in the same manner as 
those for Fig. 4. The thick dotted circle of radius r ~ 0.5Afo denotes the location of the apparent horizon. 
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FIG. 11. (a) 1 — M*(r)/M* as a function of time for r/Mo = 3 (solid curve) and 4.5 (dashed curve) for models M1616, 
M1517, and M1418. (b) The same as (a) but for models M1517 and M159183 for which the rest-mass ratio is identical but the 
total baryon rest-mass is different. 



19 




M1616 



Sf 0.6 



x 

>T 0.4 



0.2 



M1591B3 



M1418I 



I i 



M1517 



i i i i — 



1.2 



1.25 
t / P. 



1.3 



FIG. 13. Mass of apparent horizon as a function of time for models M1616 (solid curve), M1517 (dashed curve), M1418 
(long-dashed curve), and M159183 (dotted curve). 

In Fig. 13, we display the time evolution of mass of apparent horizons Mah in units of M for models M1616, 
M1517, M1418, and M159183. M AH is defined by 
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M AH = V^ (4-4) 

where S is the area of the apparent horizon. The figure indicates that Mah/M appears to approach <~ 0.85 for 
models M1616, M1517, and M159183, and ~ 0.75 for model M1418. For models of the smaller rest-mass ratios, the 
value of M An/Mo is smaller; i.e., a larger fraction of the mass element is not swallowed into a black hole. 

Since most of the mass elements fall into the black hole and radiated energy of gravitational waves is less than 
1% of M for model M1616 (see Sec. IV C), the black hole mass may be approximated by M within <~ 1% error. 
Furthermore, recall that the area of a Kerr black hole of mass M and Kerr parameter Mq (q = J/M 2 ) is written by 



S = 8nM 2 



1 + (1 - q 2 ) 1 ' 2 ] . (4.5) 



For model M1616, M ~ M and Mah ~ 0.85M . This implies a value of ~ 0.9 for the nondimensional angular 
momentum q of the formed black hole for model M1616. Since the initial value of the system, qo, is about 0.913, and 
M and J decrease by ~ 0.5% and ~ 7% by gravitational radiation (see Sec. IV), respectively, the expected final value 
of q is ~ 0.85, which agrees with the numerical value within ^5% error. (The disagreement is due to a numerical 
error associated with insufficient grid resolution. This point is confirmed by the convergence test presented in Sec. IV 
E.) The result presented here indicates that the location and the area of the apparent horizon are determined within 
~ 5% error with the current grid resolution. 



C. Gravitational waves 



In Figs. 14-16, we present the gravitational waveforms (gauge-invariant quantities) and the accumulated energy 
and angular momentum loss by gravitational radiation as a function of the retarded time (t — r o bs)/-Pt=0- (In the 
following, the retarded time is always normalized by Pt=o-) It should be noted that Rti±t jM§ = 1 implies the values 
of rh+ and rh x along the z axis due to I = m = 2 modes are w 1.85 km, where 

h + = -L Ue-^- ), h x = ^-. (4.6) 
2r A \ sin 6 J r z sin 

Irrespective of the mass and the mass ratio of binaries, the inspiral waveforms are dominant for t — r b s & Pt=0i an d 
subsequently the merger waveforms are excited. For hypermassive neutron star formation, quasiperiodic waves are 
excited because of its quasiradial and nonaxisymmetric oscillations. For the black hole formation case, the computation 
crashed soon after the formation of apparent horizon. As a result, we were not able to compute complete gravitational 
waveforms for t — r bs <; Pt=a, for which gravitational waves would be dominated by quasinormal mode ringings [41]. 
A straightforward approach to compute such gravitational waves is to develop a black hole excision technique [42], 
by which it may be possible to continue the simulation for a long time duration even after the formation of the 
black hole. An alternative approach is to extract gravitational waves from a restricted spacetime data set using the 
so-called Lazarus technique [43]. Leaving the development of the two methods for future implementations, we focus 
our discussions below mainly on the character of gravitational waves for neutron star formation than for black hole 
formation. 



1. Hypermassive neutron star formation case 

In Fig. 14, the gauge-invariant quantities for (l,m) = (2,2), (3,3), and (2,0) for models M1315 and M1414 are 
displayed. In both cases, the inspiral waveforms are dominant for t — r Q b s Ss Pt=o an d quasiperiodic waves excited by 
nonaxisymmetric quasiperiodic oscillations are emitted after the merged objects are formed for t — r Q b s <; Pt=o- It 
is found that the waveforms of the (2,2) mode for two models are similar but a slight difference can be seen in the 
quasiperiodic oscillation for t — r Q b s ^ Pt=o- For model M1315, quasiperiodic waves appear to be mainly composed 
of a single oscillation mode. On the other hand, a non-negligible modulation can be observed in the waveforms for 
model M1414. This implies that they are composed of more than two dominant modes. 

Figure 14 shows that the gauge invariant variable for the (2,0) mode does not oscillate around zero. This is due 
to the fact that gravitational waves are extracted at a finite radius and, as a result, this variable contains nonwave 
components associated with a stationary quadrupole moment. To calculate the Fourier spectrum of gravitational 
waves, we first subtracted the stationary component from i?2o and, then, performed the Fourier transformation. In 
Fig. 17, the Fourier spectra of the gauge-invariant variables for models M1315 and M1414 are displayed. Here, 
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\Rlm{f)f\r is plotted. When taking a look at this figure, the following should be also kept in mind: (i) the spectra 
presented for the (2,2) mode are not realistic for / < /qe because the spectra of inspiraling waveforms should be 
dominant in reality as |-R22(/)/| °c /~ 1//6 for / <^ /qe; (ii) the amplitude of the peaks found at / <~ 2/qe and 3/qe 
for the (2,2) mode and at /qe for the (2,0) mode is underestimated because we stopped the simulations during the 
oscillation of the formed hypermassive neutron stars to save computational time (see discussion below). 

For model M1315, a single peak is found at a frequency w 3. 2/qe in the Fourier spectrum, while for model M1414, 
two peaks of w 2.0/qe and 2.95/qe are found for the (2,2) mode. The difference in the number of the peaks reflects the 
difference of the merger process. For model M1414, the merged object constitutes a nonaxisymmetric hypermassive 
neutron star of a double core structure, which quasiradially oscillates with a large amplitude. Therefore at least two 
modes (nonaxisymmetric and quasiradial oscillation modes) are contained. The peak at / ~ 3/qe is associated with 
the nonaxisymmetric bar-mode oscillation, while that at / <~ 2/qe is produced by a modulation due to coupling 
between the nonaxisymmetric and quasiradial oscillations because the difference of their frequencies is approximately 
equal to /qe which corresponds to the frequency of the quasiradial oscillation. These two peaks in the Fourier spectra 
have been also found in Newtonian [18] and post-Newtonian simulations [19,20]. On the other hand, for model M1315, 
the merged object forms a hypermassive neutron star of an asymmetric double core structure. In this object, the 
amplitude of the quasiradial motion is not as large as that for model M1414. As a result, the peak at / ~ 2/qe 
associated with the modulation between the nonaxisymmetric and quasiradial oscillations is not as remarkable as that 
for model M1414. This feature has been also found in the post-Newtonian study [20]. 

The frequency of the peaks for the (2,2) mode around / <~ 3/qe for model M1315 is slightly larger than that 
for M1414. This results from the fact that the maximum density (or the compactness) of the formed hypermassive 
neutron stars is larger for M1315 (see Fig. 6). This fact indicates that for the smaller rest-mass ratio, the gravitational 
wave frequency associated with the nonaxisymmetric oscillation is higher for a fixed total rest-mass of the system. 

Assuming that the total mass of the system is 2.8M©, /qe is w 750 Hz for model M1414 and w 700 Hz for 
model M1315 according to Eq. (4.1). This implies that the peaks in the Fourier spectrum appear at w 1.5 and 
2.2 kHz for M1414 and at w 2.25 kHz for M1315. These frequencies will be too high to be detected by the first 
LIGO. However, these quasiperiodic gravitational waves will be interesting targets for resonant-mass detectors and/or 
specially designed advanced interferometers such as the advanced LIGO [3]. 

It should be mentioned that the peak frequencies in the post-Newtonian simulation [20] are smaller than those 
found in our study for a given neutron star mass and radius w 15 km. This may be due to the fact that, in our 
fully general relativistic simulation, the gravity is taken into account correctly and is stronger than that in the post- 
Newtonian approximation. Consequently, the formed hypermassive neutron star is more compact and hence the 
oscillation frequency higher. 

The magnitude of the quasiradial oscillation is reflected in the amplitude of gravitational waves for the (2,0) 
mode [45]. In the early phase (t — r Q b s Pt=o), this mode is dominated by a stationary quadrupole mode which 
is not associated with gravitational waves, but after a hypermassive neutron star is formed it becomes a dominant 
component. Figure 14 shows that the amplitude of the (2,0) mode for model M1414 is larger than that for M1315 by 
a factor of <~ 2. This results from the fact that the amplitude of the quasiradial oscillation for M1414 is larger than 
that for M1315. 

The frequency of gravitational waves for the (2,0) mode is within the sensitive band of kilometer-size laser inter- 
ferometers as ~ /qe ~ 0.7(2.8M Q /Mo) kHz. Although the amplitude is ~ 5% of that of the dominant (2,2) mode, 
this mode does not damp soon as indicated in Fig. 6(a). Therefore if the cycles are accumulated using a theoretical 
template, the effective amplitude may be much larger than that for one cycle and may be as large as the amplitude of 
i?22(~ Mo/r) at / ~ /qe- Unfortunately, it is difficult to exactly estimate the effective magnitude from the present 
numerical results of finite duration. However, gravitational waves of this mode may be an interesting target even for 
the first-generation gravitational wave detectors such as the first LIGO. 



22 





FIG. 15. 7?22+r and i?33+r as a function of the retarded time for models M1616 (long-dashed curves), M1517 (dashed curves), 
M1418 (solid curves), and M159183 (dotted curves). 



2. Dependence of inspiral waveforms on mass ratios 



From Fig. 15, we find that the maximum amplitude is smaller for models of the smaller rest-mass ratios. According 
to the quadrupolc formula, the maximum amplitude for a given total mass is proportional to <5m/(1 + Qm) 2 which 
is in a small range 0.246-0.25 for 0.8 < Qm < 1- This suggests that the maximum amplitude would depend weakly 
on the value of Qm- However, the ratios of the maximum amplitude for models M1517 and M1418 to that for model 
M1616 are 0.966 and 0.909, respectively. This implies that the maximum amplitude is suppressed with the decrease 
of Qm- This results from the fact that the tidal effect plays a more important role in the close binaries of the smaller 
rest-mass ratios: Since the tidal disruption sets in at a larger orbital separation for the smaller rest-mass ratios, the 
maximum amplitude should be decreased. [Similar results are also found in Fig. 16(b) (see below).] This property 
has been reported in the Newtonian and post-Newtonian studies, too [17,19,20]. According to [19,20], the suppression 
factor is proportional to Qm (for a fixed value of M„), agreeing with our results approximately 

Another difference of gravitational waveforms between models of equal-mass and unequal-mass binaries can be seen 
in the modes of odd values of m (Figs. 14 and 16). For the merger of equal-mass binary neutron stars, the amplitude 
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for those modes is zero because of the 7r-rotation symmetry. On the other hand, it is not negligible for the mergers of 
unequal-mass binaries. However, the amplitude is at most 5% of that of the (2,2) mode. 



3. Radiated energy and angular momentum 



Figure 16 shows that in the final inspiral phase (t— r Q b s Pt=o), ~ 0.3-0.5% of the initial ADM mass and ~ 6-8% of 
the initial angular momentum are carried away by gravitational radiation. (M does not change much but J decreases 
significantly.) This implies that the nondimcnsional angular momentum parameter q decreases by ~ 5-7% due to the 
gravitational radiation. Since a large fraction of baryon mass of the system is swallowed into a black hole for models 
M1616, M1517, M1418, and M159183, the ADM mass of the black hole should be ~ M within a few percents error. 
From this fact, we may expect that the final value of q is ~ 0.95(jo within a few percents error and, hence, it is in the 
range between 0.8 and 0.9. This value is approximately consistent with the value derived from the area of apparent 
horizons computed at the termination of the simulations. 

Figure 16 also indicates that the energy and angular momentum loss by gravitational radiation decrease with the 
decrease of rest-mass ratios. The reason is that tidal disruption takes place before the orbital separation becomes as 
small as the sum of radii of two stars for the merger of unequal-mass neutron stars. The orbital separation at the tidal 
disruption is larger for the smaller rest-mass ratios for a fixed value of the total rest-mass of the system. This implies 
that the maximum value of the compactness of binary orbit is smaller for the smaller rest-mass ratios, resulting in 
that the amount of gravitational radiation becomes smaller. 

For the case of hypermassive neutron star formation, the energy and the angular momentum are carried away 
gradually due to gravitational radiation emitted by the quasiperiodic nonaxisymmctric oscillations. Since the emission 
time scale is much longer than the dynamical time scale, it is impossible to follow the longterm evolution of the 
hypermassive neutron stars to the final state. If we assume that the angular momentum is dissipated by gravitational 
waves with the same rate as that at the termination of the simulation, the angular momentum will become smaller 
than O.Uo around t <; 300Pt = o- Since the hypermassive neutron stars are supported by the centrifugal force, they 
will collapse to a black hole as a result of the angular momentum dissipation within ~ 1 sec. 

In the SPH calculations [19,20,38], the quasiperiodic oscillations of the hypermassive neutron star damp in much 
shorter time than in our numerical results. If we believe their results, the lifetime of the hypermassive neutron stars 
would be much longer. The reason for the discrepancy between our and their results is unclear. However, as far 
as our simulations are concerned, there is no reason for the damping of the nonaxisymmctric oscillation in such a 
short time scale since the emission time scale of gravitational waves is much longer than one oscillation period and 
other damping processes such as dynamical angular momentum transfer are unlikely to work efficiently. We suspect 
that damping found in previous works may be due to a spurious numerical dissipation or due to an overestimation of 
gravitational radiation damping in the post-Newtonian formalism they adopted. 
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FIG. 16. The accumulated energy and angular momentum loss by gravitational radiation as a function of the retarded time 
(a) for models M1315 (solid curves) and M1414 (dashed curves), and (b) for models M1418 (solid curves), M1517 (dashed 
curves), M1616 (long-dashed curves), and M159183 (dotted curves). 
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D. Gravitational radiation reaction 



The ADM mass M and the angular momentum J computed in a finite computational domain using Eqs. (2.36) 
and (2.37) decrease with time because of the gravitational radiation. However, conservation laws (2.43) and (2.44) 
should still be satisfied. Here, we demonstrate that they arc satisfied approximately in the present simulations. 

Figure 18 shows the time evolution of M and J (solid curves) and of the quantities defined by the following equations 
(dotted curves) for models M1414 and M1517: 

M'(t) = M - AE(t), (4.7) 
J'(t) = J - AJ(t). (4.8) 

The relations M'(t) = M(t) and J'(t) = J(t) are equivalent to the conservation of the total ADM mass energy and 
angular momentum. Figure 18 indicates that relations M = M' and J = J 1 arc satisfied within <~ 1% error except for 
the phase in which the merged object collapses to a black hole and, as a result, the grid resolution becomes too poor. 

In fully general relativistic simulations, the numerical accuracy is restricted by grid resolution and by the approxi- 
mate outer boundary conditions imposed in a local wave zone. The results presented here indicate that these errors 
are suppressed within ~ 1% error in our simulations in the absence of a black hole. (In the presence of a black hole, 
the errors increased to <~ 10% and the computation crashed.) 

The conservation of the angular momentum which holds approximately in our present simulations is a necessary 
condition for studying the formation of disks and a hypermassive neutron star supported by centrifugal force, and the 
final value of q of a black hole. The results here indicate the reliability of the numerical results on the formation of 
disks and hypermassive neutron stars, and on determination of the final value of q presented in Sec. IV B and IV C. 



25 



12 3 

t / Pt-o 







0.5 1 
t / P t =o 



(a) (b) 

FIG. 18. Evolution of the total energy and angular momentum of the system calculated by Eqs. (2.36) and (2.37) (solid 
curves) and that calculated by Eqs. (4.7) and (4.8) (dotted curves) (a) for model M1414 and (b) for model M1517. 



E. Calibrations 




FIG. 19. Time evolution of the maximum density, the central values of a and (ft, the averaged violation of the Hamiltonian 
constraint (H), M, and J for models M1616-3 (dotted curves), M1616-4 (solid curves), and M1616-5 (long-dashed curves). In 
this figure, effects with regard to the grid resolution are clarified. 
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FIG. 20. The same as Fig. 19, but for models M1616 (solid curves), M1616-2 (dashed curves), and M1616-3 (long-dashed 
curves). In this figure, effects with regard to the location of the outer boundaries are clarified. 
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(a) (b) 
FIG. 21. (a) Gravitational waveforms (R22+), and accumulated energy and angular momentum of gravitational radiation for 
models M1616 (solid curves), M1616-2 (dashed curves), M1616-3 (long-dashed curves), and M1616-4 (dotted curves), (b) The 
same as Fig. 18(b), but for models M1616 (solid and dotted curves) and M1616-3 (dashed and dotted-dashed curves). The 
solid and dashed curves denote M/Mq and J/ Jo, and the dotted and dotted-dashed curves 1 — AE/Mo and 1 — AJ/Mq. 

Convergence tests were performed employing models M1616 and M1414. The test simulations were done for 
additional five models as listed in Table II. To investigate effects of the location of the outer boundaries at which 
approximate boundary conditions were imposed, the values of L were changed for three levels as L/X — 0.533 (M1616), 
0.425 (M1616-2), and 0.263 (M1616-3) fixing the grid spacing. To see effects with regard to the grid resolution, we also 
performed two additional simulations for models M1616-4 and M1616-5 in which the location of the outer boundaries 
was the same as that for M1616-3, but the grid spacings were about 5/6 and 5/4 times, respectively, that for M1616- 
3. A simulation for model M1414-2 was performed to clarify weak dependence of gravitational waveforms from 
quasiperiodic oscillations of a hypermassive neutron star on the value of L. 



1. Convergence test with regard to grid resolution 

Figure 19 shows the evolution of the maximum density, the central values of a and <j>, the averaged violation of 
the Hamiltonian constraint H [computed by Eq. (2.46)], M [computed by Eq. (2.36)], and J [computed by Eq. 
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(2.37)] for models M1616-3 (dotted curves), M1616-4 (solid curves), and M1616-5 (long-dashed curves). This figure 
indicates the dependence of the numerical results on the grid resolution for a fixed value of L. It is found that the 
convergence of H is at first order. A likely reason is as follows: Since the vacuum is not allowed in our hydrodynamic 
implementation, we have to add an atmosphere of small density outside neutron stars. In the present work, the 
density of the atmosphere is <~ 10~ 7 in units of k = 1. As a result, a very steep density gradient appears at the stellar 
surface. In such a region, the transport term of the hydrodynamic equations is computed with first-order accuracy in 
space. This effect seems to be non-negligible in determining the global order of the accuracy. 

The angular momentum is dissipated and transported unphysically by numerical effects. For the larger grid spacing, 
the dissipation rate is larger and, as a result, the duration of the inspiral phase becomes shorter. Even in the case 
of the best resolution (M1616-4), the angular momentum appears to be dissipated by ~ 1%. This effect may be the 
main source for the discrepancy between J and J' [see Fig. 18(a)] in the late phase t — r obs <; P t=0 for model M1414. 



2. Convergence test with regard to L 

In Fig. 20, we show the same figure as that of Fig. 19 but for models M1616 (solid curves), M1616-2 (dashed 
curves), and M1616-3 (long-dashed curves) to make a comparison among the numerical results with the different 
values of L and a fixed grid spacing. It is found that (i) for a smaller value of L, the merged object collapses earlier, 
(ii) H depends very weakly on the value of L, and (iii) the results for models M1616 and M1616-2 are almost identical. 

The reason for (i) is that the magnitude of the radiation reaction is overestimated with small values of L. To explain 
this effect, gravitational waveforms, the radiated energy, and the radiated angular momentum are shown in Fig. 21, 
which indicate that the numerical results for models M1616 and M1616-2 are approximately identical. This implies 
that with L J> 0.5Ao, a convergent result may be achieved. On the other hand, with the smaller value of L < 0.5Ao, the 
amplitude of gravitational waves, the radiated energy, and the radiated angular momentum are overestimated. The 
radiated energy and angular momentum for model M1616-3 are about twice as large as those for M1616. As a result, 
the orbital separation for M1616-3 decreases more rapidly than that for M1616. Moreover, the radiation reaction is 
not accurately computed for model M1616-3, so that the conservation of the angular momentum (J + AJ = Jo) is 
largely violated [see Fig. 21(b)]. A number of numerical simulations for the binary merger in full general relativity 
have been recently performed with L < 0.5Ao [7,11,12]. Figure 21(b) warns that the gravitational waveforms and the 
merger process in such numerical simulations are not very reliable. 




FIG. 22. Gravitational waveforms (R22+) and accumulated energy and angular momentum of gravitational radiation for 
models M1414 (solid curves) and M1414-2 (dashed curves). 

Figure 22 is the same figure as that of Fig. 21(a) but for models M1414 (solid curves) and M1414-2 (dashed curves), 
for which the outer boundaries are located at L = 0.510A and 0.252A , respectively. For t— r obs <J P t =o, the amplitude 
of gravitational waves, the radiated energy, and the radiated angular momentum are overestimated for the smaller 
value of L. Since the angular momentum is dissipated more rapidly from the system, the inspiral phase is shorter 
and the merger sets in earlier for model M1414-2. This results in a phase difference between gravitational waves of 
M1414 and M1414-2. However, for t — r Q bs ^ Pt=o, the amplitude of gravitational waves, the energy luminosity, and 
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the angular momentum flux are approximately in agreement between two results (besides the slight disagreement in 
the wave phase). This figure shows that quasiperiodic waves emitted from oscillating hypermassive neutron stars is 
calculated accurately with our choice of L, since its wavelength is short enough to compute these quantities even for 
the smaller value of L. 



V. SUMMARY 

We performed fully general relativistic simulations for the merger of binary neutron stars focusing particularly on 
the unequal-mass case. The following is a summary of the scientific results obtained in this paper: 

• If the total rest-mass of the system is more than 1.7 times of the maximum allowed rest- mass of spherical neutron 
stars a black hole is formed for the T-law equation of state with n = 1. The nondimcnsional angular momentum 
parameter of the formed Kerr black hole is likely to be in the range between 0.8 and 0.9. 

• Disk mass around a black hole formed after the merger increases with the decrease of rest-mass ratios for a fixed 
value of the total baryon rest-mass of binary neutron stars. It is found that for the rest-mass ratio ~ 0.85, the 
disk mass may be several percents of the total mass of the system if two neutron stars are not very compact. 

• Disk mass around a black hole formed after the merger decreases with the increase of the compactness of the 
system for a fixed value of the rest-mass ratio. 

• Shape of the hypermassive neutron stars formed after the merger depends on the rest-mass ratio of binaries. 
For the merger of equal-mass neutron stars, a hypermassive neutron star of a double core is formed. On the 
other hand, for the merger of unequal-mass neutron stars, an asymmetric double core structure is the outcome. 

• In the hypermassive neutron stars formed after the merger, both nonaxisymmetric and quasiradial oscillations 
are excited. These oscillations induce gravitational radiation. 

• For the case of hypermassive neutron star formation, the characteristic frequency of gravitational waves associ- 
ated with nonaxisymmetric oscillations is ~ 3/qe, which is w 2.2 kHz assuming that Mq w 2.8Mq. This value 
is slightly higher than that found in the post-Newtonian simulation [20]. This is likely due to the fact that the 
formed hypermassive neutron star is more compact in our simulation in which general relativistic effects are 
fully taken into account. 

• The frequency of the peak in the gravitational wave spectrum associated with the nonaxisymmetric oscillation 
is higher for the mergers of the smaller rest-mass ratio with a given total rest-mass. This reflects the fact that 
the formed hypermassive neutron star is more compact for mergers of the smaller rest-mass ratios. 

• The amplitude of quasiradial oscillations for hypermassive neutron stars is larger for the merger of equal-mass 
neutron stars. This is reflected in the amplitude of gravitational waves for R20 as well as the magnitude of the 
peak at ~ 2/ QE of R 2 2 ■ 

• The characteristic frequency of gravitational waves associated with a quasiradial oscillation is ~ /qe- The 
oscillation does not damp quickly. Thus if the cycle of gravitational waves could be accumulated using a 
theoretical template, the effective amplitude may be as large as that of the dominant quadrupole component. 

The simulations were performed using a new implementation. As a result, the accuracy of the numerical results is 
significantly improved. In particular, we emphasize that the gravitational radiation reaction is taken into account with 
a good accuracy in the new implementation. We now consider that fundamental parts of the numerical implementation 
such as those for Einstein's evolution equation, general relativistic hydrodynamic equations, gauge conditions, and 
apparent horizon finder are established well for simulating spacetimes of no black hole and for early growth of formed 
black holes. 

However, there are still technical issues to be solved. The following is a list of them: 

• The black hole forming region does not have good resolution in our current computation. Consequently, com- 
putation crashed soon after formation of the apparent horizon. Obviously, it is necessary to improve the grid 
resolution around the black hole forming region for longer time simulations. Since we have to prepare a large 
computational domain with L which is at least half of the wavelength of gravitational waves, using restricted 
computational speed and memory, it is desirable to develop numerical techniques such as the mesh refinement 
techniques [44] to overcome this problem. 
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• Gravitational waveforms are incompletely computed in the case of black hole formation, since the computations 
crash soon after the formation of the black holes. A straightforward approach to compute such gravitational 
waves is to develop a black hole excision technique [42] by which we might be able to continue the simulation for 
a long time duration even after formation of the black holes. An alternative approach is to extract gravitational 
waves from a restricted spacetime data set using the so-called Lazarus technique [43] . Developing either of two 
technologies is an issue for the future. 

• Up until this time, we have performed simulations using T-law equations of state and neglecting micro physical 
effects. To produce more physical and realistic outputs by numerical simulation, it is necessary to take into 
account sophisticated microphysics as done in Newtonian simulations [46]. 

• It is desirable to improve the implementation for providing the initial conditions. In simulations performed to 
this time, we have used quasiequilibrium states of a conformally flat three-metric as the initial conditions for 
simplicity. The conformal flatness approximation becomes a source of a certain systematic error when attempting 
to obtain realistic quasiequilibrium states, since the nonconformal part of the three-metric is in general nonzero 
[47]. As a result, this approximation introduces a systematic error on the initial conditions and subsequent 
merger simulation. Since the magnitude of the ignored terms in the conformal flatness approximation seems to 
be small, it is unlikely that this effect significantly changes the results obtained in this paper. However, this 
conclusion is not entirely certain. To rule out the possibility, it is necessary to perform simulations using 
quasiequilibrium states of generic three geometries as initial conditions. A few formulations in which the 
conformal flatness is not assumed have been proposed recently [48] . 
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Model 


{M/R) x 


Pmax 


Qm 


M t 


M 


<7o 


Pt=o 


Co 


Q, 


L/Ao 


Product 


M disk /M» 


M1414 


0.14 vs 0.14 


0.118, 0.118 


1.00 


0.292 


0.269 


0.951 


193 


0.102 


1.62 


0.510 


NS 




M1315 


0.13 vs 0.15 


0.104, 0.134 


0.901 


0.292 


0.269 


0.961 


206 


0.0976 


1.62 


0.479 


NS 




M1616 


0.16 vs 0.16 


0.151, 0.151 


1.00 


0.320 


0.292 


0.914 


158 


0.116 


1.78 


0.533 


BH 


< 0.2% 


M1517 


0.15 vs 0.17 


0.133, 0.171 


0.925 


0.319 


0.291 


0.923 


169 


0.111 


1.77 


0.507 


BH 


£2% 


M1418 


0.14 vs 0.18 


0.118, 0.195 


0.855 


0.317 


0.290 


0.933 


182 


0.106 


1.76 


0.467 


BH 


£4% 


M159183 


0.159 vs 0.183 


0.149, 0.203 


0.925 


0.332 


0.301 


0.908 


156 


0.118 


1.84 


0.498 


BH 


< 1% 



TABLE I. A list of several quantities for quasiequilibria of irrotational binary neutron stars with n = 1. The compactness 
of each star in isolation (M/R)^, the maximum density for each star, the baryon rest-mass ratio Qm = M^/M*!, the total 
baryon rest-mass, the ADM mass at t = (Mo), go = Jo/Mq, Pt=o = Pt=o/Mo, the orbital compactness [Co = (Mofi) 2 ' 3 ], 
the ratio of the total baryon rest-mass to the maximum allowed mass for a spherical star (Q, = M»/M| P m ax ), gravitational 
wavelength in units of L in the maximum grid number, and the products we found when we stopped simulations, In the last 
column, the estimated ratio of the disk rest-mass located for r > 3Mo at the termination of the simulation to the total rest-mass 
for the black hole formation case is listed. All quantities are normalized by n appropriately to be dimensionless: The mass, the 
radius, and the density can be rescaled to desirable values by appropriately choosing k. Here, M, sp J,„ denotes the maximum 
allowed mass of a spherical star (M* P m ax ~ 0.180 at p max ~ 0.32 for n — 1 and k = 1). BH and NS denote "black hole" and 
"neutron star" . 



Model 


A/Mo 


Grid size 


L/\ 


L/M 


M1616 


0.134 


(633,633,317) 


0.533 


42.2 


M1616-2 


0.134 


(505,505,252) 


0.425 


33.7 


M1616-3 


0.134 


(313,313,157) 


0.263 


20.8 


M1616-4 


0.111 


(377,377,189) 


0.263 


20.8 


M1616-5 


0.169 


(249,249,125) 


0.263 


20.8 


M1414 


0.156 


(633,633,317) 


0.510 


49.3 


M1414-2 


0.156 


(313,313,157) 


0.252 


24.3 



TABLE II. Computational setting for test simulations. 
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